/*
 * AqueousActivity.cpp
 *
 *  Created on: 23 May 2011
 *      Author: allan
 */

#include "AqueousActivity.h"

// C++ includes
#include <algorithm>
using namespace std;

// GeoReact includes
#include "../AqueousSolution.h"
#include "Aqueous/HkfModel.h"
#include "Aqueous/SetschenowModel.h"
#include "Aqueous/DrummondCO2Model.h"

const AqueousActivityModel IdealModel()
{
	auto waterMoleFraction = [](const AqueousSolution& sol, double T, double P, const VectorXd& n) 
			{ return n[sol.GetWaterIndex()] / n.sum(); };
	
	return waterMoleFraction;
}

const AqueousActivityModel HkfIonicModel(const string& ion, const vector<string>& ions)
{
	vector<double> effRadii;
	
	for(unsigned i = 0; i < ions.size(); ++i) 
		effRadii.push_back(EffectiveElectrostaticRadius(ions[i]));
	
	const Index indexIon = std::find(ions.begin(), ions.end(), ion) - ions.begin();
	
	return std::bind(HkfIonicActivityCoefficient, _1, _2, _3, _4, indexIon, effRadii);
}

const AqueousActivityModel HkfWaterModel(const vector<string>& ions)
{
	vector<double> effRadii;
	
	for(unsigned i = 0; i < ions.size(); ++i) 
		effRadii.push_back(EffectiveElectrostaticRadius(ions[i]));
	
	return std::bind(HkfWaterActivityCoefficient, _1, _2, _3, _4, effRadii);
}

const AqueousActivityModel SetschenowModel()
{
	return SetschenowModel(0.0);
}

const AqueousActivityModel SetschenowModel(double b)
{
	return SetschenowModel({b}, {25.0}); // Assume the value b is given at 25.0 degC
}

const AqueousActivityModel SetschenowModel(const vector<double>& bi, const vector<double>& Ti)
{
	return std::bind(SetschenowActivityCoefficient, _1, _2, _3, _4, bi, Ti);
}

const AqueousActivityModel DrummondCO2Model()
{
	return AqueousActivityModel(DrummondCO2ActivityCoefficient);
}
